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ABSTRACT 

We analyse parallel N-body simulations of three Cold Dark Matter (CDM) universes 
to study the abundance and clustering of galaxy clusters. The simulation boxes are 
500/i _1 Mpc on a side and cover a volume comparable to that of the forthcoming Sloan 
Digital Sky Survey. The use of a treecode algorithm and 47 million particles allows us 
at the same time to achieve high mass and force resolution. We are thus able to make 
robust measurements of cluster properties with good number statistics up to a redshift 
larger than unity. We extract halos using two independent, public domain group find- 
CTN ' ers designed to identify virialised objects - 'Friends-of-Friends' (Davis et al. 1985) and 

'HOP' (Eisenstein & Hut 1998) - and find consistent results. The correlation function 
of clusters as a function of mass in the simulations is in very good agreement with 
I a simple analytic prescription based upon a Lagrangian biasing scheme developed by 

Mo & White (1996) and the Press-Schechter (PS) formalism for the mass function. 
The correlation length of clusters as a function of their number density, the Rq-D c 
relation, is in good agreement with the APM Cluster Survey in our open CDM model. 
The critical density CDM model (SCDM) shows much smaller correlation lengths than 
are observed. We also find that the correlation length does not grow as rapidly with 
cluster separation in any of the simulations as suggested by the analysis of very rich 
Abell clusters. Our SCDM simulation shows a robust deviation in the shape and evo- 
lution of the mass function when compared with that predicted by the PS formalism. 
Critical models with a low as normalization or small shape parameter T have an excess 
of massive clusters compared with the PS prediction. When cluster normalized, the 
SCDM universe at z = 1 contains 10 times more clusters with temperatures greater 
than 7keV, compared with the Press & Schechter prediction. The agreement between 
the analytic and N-body mass functions can be improved, for clusters hotter than 3 
keV in the critical density SCDM model, if the value of 5 C (the extrapolated linear 
theory threshold for collapse) is revised to be 5 c (z) = 1.685 [(0.7/<78)(l + ' (as 
is the rms density fluctuation in spheres of radius 8/i _1 Mpc). Our best estimate for 
the amplitude of fluctuations inferred from the local cluster abundance for the SCDM 
model is as — 0.5 ± 0.04. However, the discrepancy between the temperature func- 
tion predicted in a critical density universe and that observed at z = 0.33 (Henry et 
al. 1998) is reduced by a modest amount using the modified Press-Schechter scheme. 
The discrepancy is still large enough to rule out Qq = 1, unless there are significant 
differences in the relation between mass and temperature for clusters at high and low 
redshift. 

Key words: cosmology- clusters- general- large scale structure of the universe. 



in 
> 



© 1998 RAS 



2 F.Governato et al. 



1 INTRODUCTION 

Clusters of galaxies, by virtue of being both relatively rare 
objects and the largest gravitationally bound systems in the 
Universe, provide stringent constraints on theories of struc- 
ture formation. The two cluster properties that are most 
commonly discussed in this context are the abundance and 
the spatial clustering. The model predictions depend sen- 
sitively on the cosmology and on the value of ag, the rms 
density fluctuations on the scale of 8 hT 1 Mpc. (Here and 
throughout this paper, h is the present-day Hubble constant 
in units of 100 km/s/Mpc.) Comparisons between observa- 
tions and model predictions have been used to place con- 
straints on cosmological parameters (Strauss et al 1995, Eke, 
Cole & Frenk 1996, Viana & Liddle 1996; Mo, Jing & White 
1996, Borgani et al 1997, De Theije, Van Kampen & Sli- 
jkhuis 1998, Postman 1998 ). 

It has long been known that clusters of galaxies are 
much more strongly clustered than galaxies (see, for exam- 
ple, Hauser & Peebles 1973 and review by Bahcall 1988). The 
two-point correlation function for the clusters is roughly a 
power law: £ cc (r) = (r/_R )" 1,8 - Bahcall & West (1992) ar- 
gue that the correlation length, _Ro, obeys the scaling rela- 
tion 

Ro&QAD c , 20/T 1 Mpc < D c < 100ft -1 Mpc, (1) 

— 1/3 

where D c = n c is the mean intercluster separation and 
n c , is the mean space density of clusters. The combined set 
of results based on the analysis of the spatial clustering of 
an X-ray flux-limited sample of clusters (Lahav et al. 1989; 
Romer et al. 1994, Abadi, Lambas & Muriel 1998), of clus- 
ters containing cD galaxies (West & Van den Bergh 1991), 
of richness class R>0, R>1, R>2 Abell clusters (Pea- 
cock & West 1992; Postman, Huchra & Geller 1992), and 
of the cluster samples extracted from the APM Galaxy Sur- 
vey (Dalton et al. 1992) and Edinburgh-Durham Southern 
Galaxy Catalogue (Nichol et al. 1992) all give results that 
are roughly consistent with the above scaling relation. 

However, on scales greater than D c ~ 40ft" 1 Mpc, the 
evidence in favour of the scaling relation hinges just on 
the analyses of the R > 1 and R > 2 Abell cluster sam- 
ples, which give _Ro ~ 21ft -1 Mpc for D c ~ 55/t -1 Mpc and 
Ro ~ 45/i _1 Mpc for D c fa 94/i _1 Mpc, respectively. Several 
authors (e.g. Sutherland 1988; Dekel et al. 1989; Sutherland 
& Efstathiou 1991) have suggested that these correlation 
lengths have been biased upward by the inhomogeneities 
and projection effects in the Abell catalogue. However, this 
suggestion has been rejected by, for example, Jing, Plionis 
& Valdarnini (1992) and Peacock & West (1992). More re- 
cently, Croft et al. (1997) have analyzed the correlation 
properties of a sample of "rich" APM clusters and find that 
the cluster correlation length saturates at Ro « 15/i _1 Mpc 
(Ro ~ 20h~ 1 Mpc if the analysis is done in redshift-space — 
see Croft et al. 1997) for D c > 40/i _1 Mpc. The controversy 
regarding the correlation length of rich clusters: i.e. if the 
i?o vs D c flattens at large scales is as of yet still unresolved. 

In an effort to resolve this issue, several authors 
(e.g. Bahcall & Cen 1992; Watanabe et al. 1994; Croft & 
Efstathiou 1994, 1997; Walter & Klypin 1996, Eke et al. 
1996) have turned to large numerical simulations. Bahcall 
& Cen (1992) investigated the cluster correlation properties 
in large N-body simulations of the standard CDM model 



(SCDM) and two low-f^o models ( Qo is the density pa- 
rameter), one spatially flat and one open. They claim to 
find a linear relation between Ro and D c over the range 
30/i _1 Mpc < D c < 95ft" 1 Mpc in all the models but that 
only in the low-density models is the Ro-D c relation steep 
enough to be consistent with the suggested scaling relation 
(0). More recent works (Croft & Efstathiou 1994, Watan- 
abe et al. 1994) have confirmed that the SCDM model is 
incompatible with the observed degree of clustering on all 
scales and for all normalizations. However, no general agree- 
ment was reached on the clustering strength at large scales 
for the other models investigated. 

In summary, apart from the general agreement that the 
SCDM model fails to account for the observed cluster cor- 
relations, results obtained from the numerical studies, due 
to lack of consistency, have been singularly unhelpful in re- 
solving the cluster correlation controversy. 

If cluster correlations are going to be used to constrain 
models of structure formation and place limits on the values 
of the fundamental cosmological parameters, it is important 
to understand why these numerical studies give such dis- 
crepant results. This is a necessary step before a meaningful 
comparison between theoretical (numerical) predictions and 
observations is possible. There are several factors that can 
affect numerical results and cause the discrepancy described 
above. Among these are differences in the mass and force 
resolution of the simulations as well as the overall volume of 
the simulations. Rich clusters tend to be rare objects and, 
therefore, simulation studies of the properties of such ob- 
jects must necessarily span large cosmological volumes. Of- 
ten, computational limitations require that such simulation 
studies compromise on the resolution (mass and/or force). 
However, this can have serious effects on the results. Watan- 
abe et al. (1994), have shown that degrading the mass reso- 
lution tends to bias the correlation lengths downward. Con- 
sequently, there is a definite need for analysis of a sample 
of simulated clusters extracted from a simulation with high 
mass and force resolution, large number of time steps, and 
covering a sufficiently large cosmological volume. 

In addition to differences in resolution and size, there 
is the issue of how to identify clusters in the simulations. 
Bahcall & Cen (1992), Watanabe et al. (1994) and Croft et 
al. (1994; 1997) used different algorithms to identify clus- 
ters in their simulations. Using a O = 1 SCDM model Eke 
et al. (1996b) explored the possibility that different algo- 
rithms could indeed give different results. They identified 
and ranked the clusters in their simulations in various dif- 
ferent ways, and found that each algorithm/selection criteria 
imprints its own particular set of biases on the cluster sam- 
ple; for a fixed value of D c , the clustering length can vary 
up to a factor of ~ 1.5. 

In this paper, we report on our analysis of cluster cor- 
relations in simulations of both critical density (Q — 1) and 
open (fio = 0.3 and fio = 0.4) CDM cosmogonies and use 
our results to explore the questions raised above. As de- 
scribed in the next section, both the force and mass reso- 
lution of our simulations are better than those of previous 
studies. Moreover, our simulated volumes are comparable 
with the Sloan Digital Sky Survey (Loveday 1998) and are 
larger than the 2dF survey (Colless & Boyle 1998). 

We investigate the present-day abundances and tem- 
poral evolution of the abundances in the three CDM mod- 
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els. Specifically, we are interested in testing the validity of 
the widely-used analytic Press-Schechter expression for the 
cluster mass function. The combination of the present-day 
abundance of clusters and the rate at which the abundance 
evolves as a function of time place strong constraints on fio 
and erg. (White, Efstathiou & Frenk 1993; Viana & Liddle 
1996; Eke, Cole & Frenk 1996). Since real clusters of galaxies 
are the product of non-linear gravitational and gas dynami- 
cal processes, the most direct way of constraining the range 
of fio and erg is to carry out large-scale numerical simula- 
tions of different models that have the necessary dynami- 
cal range and include (poorly known) gas-stellar physics, 
then "observe" the resulting model universe and compare 
the simulated observations with the real ones. Computa- 
tionally, this route is prohibitively expensive at present. A 
more economical approach involves using the analytic Press- 
Schechter (PS) formalism (Press & Schechter 1974; Bond et 
al. 1991) to compute the cluster mass function, map the 
mass function into an abundance distribution as a function 
of the observable parameter, and then compare the latter to 
observations in order to determine the appropriate values of 
fio and erg. Setting aside the uncertainties in the correspon- 
dence between mass and an observable quantity, the validity 
of the analytic approach rests entirely on the assumption 
that the PS formalism yields an accurate description of the 
cluster mass function. The analytic expression for the cluster 
mass function has been extensively tested against numerical 
simulations in the past (see, for example, Carlberg & Couch- 
man 1989; Lacey & Cole 1994; Klypin et al. 1994, Cole, 
Weinberg, Frenk & Ratra 1997, Cen 1998) and most studies 
have found a good agreement between the analytical and 
the numerical results. However, there have also been some 
interesting claims to the contrary. Gross et al. (1998), for 
example, have drawn attention to a discrepancy between the 
PS predictions and numerical results at small masses, and 
Bertschinger & Jain (1994) claim that the PS mass function 
systematically underestimates the number density of high 
mass halos. Estimates of parameters such as fio and erg are 
usually derived from fitting the analytic cluster mass func- 
tion to the observed distribution. If the PS mass function is 
indeed failing at the high mass end and this failure is not 
taken in account, it can affect the determinations of fio and 

In this paper our aim is to determine the halo mass 
function on group and cluster scales in our set of simula- 
tions, and use these to assess the reliability of the analytic 
PS mass function. Each of our volumes contains several hun- 
dred 'Coma-like' clusters at the present time. This, in con- 
junction with our high mass and force resolution, allows us 
to map out the cluster mass function to high precision out 
to z ~ 1, i.e. over a larger redshift range than previously 
possible. 

The lay-out of the paper is as follows: In §2, we discuss 
our numerical simulations and the procedure for construct- 
ing cluster catalogs. We use these catalogs to extract the 
cluster mass function and to study their spatial correlation 
properties. In §3, we present the results of our correlation 
analyses and in §4, we discuss the cluster mass function. In 
both sections we compare our numerical results to analyt- 
ical approximations. Adopting a simple mapping between 
mass and X-ray temperature, we transform our numerical 
mass function into a temperature function and highlight the 



main differences between this temperature function and the 
one based on the standard Press-Schechter mass function. 
Finally, we summarize our results and briefly discuss their 
relevance for future cosmological tests in §5. 



2 NUMERICAL SIMULATIONS AND 
CLUSTER SELECTION 

We have simulated structure formation within a periodic 
cube of comoving length L = 500ft" 1 Mpc for two "fidu- 
cial" cosmological models: A critical density cold dark mat- 
ter models (fio = 1, h = 0.5 with erg = 1.0 at z = ( here- 
after we refer to the erg = 0.7 output as SCDM07 and to 
the erg = 1 as SCDM10, respectively. Of course each output 
of the SCDM run can be rescaled to a different z changing 
the present day erg normalization) and an open (fio = 0.3, 
h — 0.75, erg = 1.0 at z = — hereafter referred to as 
03CDM) cold dark matter model. The z = 0.58 output 
of 03CDM simulation can, with appropriate rescaling, be 
identified as an fi = 0.4, h = 0.65, erg = 0.79 CDM simula- 
tion (hereafter referred to as 04CDM) of comoving length 
L = 433.3 h" 1 Mpc. The same set of simulations was used 
by Szapudi et al. (1998) to study the higher order corre- 
lation properties of galaxies. The initial conditions were set 
using the Bardeen et al (1986) transfer function for CDM. 
The simulations were computed using PKDGRAV, a par- 
allel treecode that allows for periodic boundary conditions 
and individual time steps (Stadel & Quinn, in preparation). 
These are among best studied cosmological models (Davis 
et al. 1985, Jenkins et al. 1997), our choices for the normal- 
ization (erg) of the open models correspond roughly to those 
inferred from the present-day cluster abundance (see, e.g., 
Eke et al. 1996 and references therein), while we analyzed 
the SCDM simulation data on a erg range that goes from erg 
= 1 (roughly COBE normalized) to erg = 0.35 (correspond- 
ing to z = 1.85 for a erg = 1 at z = model, and to z = 0.43 
for a cluster normalized SCDM universe with erg = 0.5 at z 
= 0) . A cubic spline force softening of 50/i _1 kpc (43/i _1 
kpc for 04CDM) was used so that the overall structure of 
clusters could be resolved. Accurate forces were maintained 
by using a cell opening angle of 8 < 0.8 (or better at high 
z) and expanding the potentials of cells to hexadecapole or- 
der. Timesteps were constrained to At < 0.3a/ e/o, where 
e is the softening length and a is the magnitude of the ac- 
celeration of a given particle. See Quinn et al. (1997) for a 
discussion and tests of this timestep criterion. In each run 
47 million particles were used, arranged on a 360 3 grid. Each 
run took several hundred hours on 256 nodes of a Cray T3E 
supercomputer, and about a thousand timesteps. The parti- 
cle mass is 7.4xlO 11 r7fi o /i- 1 M , where rj = 1 for SCDM and 
03CDM models and r\ = 0.65 for 04CDM. Simulations were 
started at z = 49. The extremely large volumes simulated, 
coupled with a reasonable mass resolution and the very good 
force resolution made possible by the use of a treecode, allow 
us to study in detail the evolution of structures ranging in 
size from groups of galaxies, made up of several tens of par- 
ticles each, to very rich clusters that contain a few thousand 
particles. In our analyses, we only consider halos consist- 
ing of 64 particles or more. This is a stricter constraint than 
used in most previous work, and was imposed to ensure that 
our results were not influenced by small-number effects. Fi- 
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nally, we verified that in these simulations both the initial 
and present-day power spectrum were in close agreement 
with theoretical expectations ( see Peacock & Dodds 1996) . 

2.1 Cluster Identification and Selection 

Theoretical treatments generally define virialized halos at 
a given epoch as structures with a mean density averaged 
over a sphere of ~ 200 times the critical density at that 
epoch (see, for example, Lacey & Cole 1994 and references 
therein). The mass contained in the sphere is taken to be 
the mass of the halo and the radius of the sphere is usually 
identified as the virial radius of the halo. 

In numerical simulations, halos are identified using a 
variety of schemes. Of these, we have chosen to use two that 
are available in the public domain: FOfQ (Davis et al. 1985), 

and HOP0 (Eisenstein & Hut 1998). These schemes are dis- 
cussed in the next subsection. Other halo finders that are of- 
ten used in literature to find virialized halos are DENMAX 
(Gelb & Bertschinger 1994), the "spherical overdensity al- 
gorithm " or SO, that finds spherically averaged halos above 
a given overdensity (Lacey & Cole 1994) and the scheme re- 
cently developed by Gross et al (1998). The algorithms that 
we opted to use are those in the public domain and hence, 
in common use. We felt that it was important to ascertain 
the extent to which these schemes may bias our results. 

2.2 Friend of Friends: FOF 

The FOF algorithm (Davis et al. 1985) is one of the most 
widely used. It is based on a nearest neighbor search. The 
main advantages of this algorithm are its simplicity and the 
lack of assumptions about the shape of halos. In this scheme, 
all particle pairs separated by less than b times the mean in- 
terparticle separation are linked together. Sets of mutually 
linked particles form groups that are then identified as dark 
matter halos. In the present study, we adopted the linking 
length that Lacey & Cole (1994) arrived at to identify viri- 
alized halos with mean densities of ~ 200 times the critical 
density at the epoch under consideration. The linking length 
is 0.2f2(z) -1//3 times the mean comoving interparticle sep- 
aration. Moreover in the low Q models, the scaling of the 
linking length as a function of redshift was further modi- 
fied as the mean halo density associated with virialization 
is a function of redshift (see for example, Kitayama & Suto 
1996). The resulting halos also have a mass-radius relation 
that agrees reasonably well with the theoretical relation for 
virialized halos (see Lacey & Cole 1994; also, Eke, Cole, 
Frenk & Navarro 1996). The objects identified by the FOF 
algorithm are the kinds of objects that the PS formalism 
refers to (apart from the lack of spherical symmetry) and 
therefore, we should be able to make a meaningful compari- 
son between the distribution of halos in the simulations and 
the PS distribution. Several authors have reported the ten- 
dency of FOF to link together close binary systems of similar 
mass especially if the two happen to be loosely connected 
by a bridge of particles. This pathology can, in specific cases 

* http:/ /www-hpcc. astro. washington.edu/tools/FOF/ 
t http:/ /www. sns.ias.edu/eisenste/hop/hop. html 



(see for example Governato et al. 1997), give rise to biased 
results. We have verified that our results are largely unaf- 
fected by this problem. 

We note that in their study, Lacey & Cole (1994) com- 
pared the properties of the halo population defined by FOF 
and with those of a sample generated using the SO algo- 
rithm. They found that at least over the mass range that 
they could probe using their simulations, the two algorithms 
gave very similar results. 

2.3 HOP 

HOP is a recently introduced algorithm (Eisenstein & Hut 
1998) based on an hybrid approach. The local density field is 
first obtained by smoothing the density field with an SPH- 
like kernel using the n nearest neighbours (we used 16) . The 
particles above a given threshold are linked with their high- 
est density neighbors until, after several "hops", they are 
connected to the one particle with the highest density within 
the region above the threshold. All particles linked to the lo- 
cal density maximum are identified as a group. Like FOF, 
HOP is well suited to identifing virialized structures once the 
density threshold is specified to be the local density at the 
virial radius. Eisenstein & Hut (1998) claim good agreement 
with the FOF method at masses above the smoothing scale. 
However, HOP can be tuned to separate binary halos — 
binary systems loosely connected by one-dimensional parti- 
cle bridges — thereby avoiding the (rare) FOF pathology. 

2.4 FOF vs. HOP 

In Figure [j] we show the mass functions obtained by apply- 
ing the two halo finders described to the 03CDM run. The 
results for SCDM runs are qualitatively similar. The FOF 
and HOP mass functions agree quite well over the entire 
mass range probed, with most massive HOP clusters show- 
ing a systematic offset of about 7% toward larger masses. 
This offset can be easily adjusted changing b or the density 
threshold for HOP. However, as discussed in the above para- 
graph the parameters used are the most physically meaning- 
ful and the small offset is a measure of the kind of biases 
that you get using different halo finders. 

We have found this general agreement to hold for dif- 
ferent models and at all redshifts. This result, coupled with 
results of Eke et aL's (1996) comparison of the FOF and 
SO algorithms, strongly indicates that regardless of the ac- 
tual details of the scheme used to identify the halos, if the 
resulting halos are independent virialized entities then the 
statistical properties of the halo populations will be very 
similar. Figure |l] also shows the PS prediction as a compar- 
ison. We defer the comparisons of the theoretical curve to 
the numerical results in §4. 



3 THE TWO-POINT CLUSTER 
CORRELATION FUNCTION 

The output of our cosmological simulations was processed 
using the two halo identification algorithms (HOP & FOF) 
described in the preceding section. We ordered the lists ac- 
cording to halo mass and then generated cluster catalogs by 
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M > 4.6xl0 13 h- 1 M 

• SCDM07-F (I3h- 1 ) 
OSCDM07-H (13h-') 
A 03CDM -F (21b- 1 ) 
□ 03CDM -H (20b.- 1 ) 



D c = 15 h- 1 Mpc 

• SCDM07-F 
O SCDM07-H 
A 03CDM -F 
□ 03CDM -H 





M > 2.5xl0 14 h- 1 M 

• SCDM07-F (28h-'; 
OSCDM07-H (26h-') 
A 03CDM -F (47h- r 
□ 03CDM -H (42h- ] 



• SCDM07-F 
OSCDM07-H 
A 03CDM -F 
□ 03CDM -H 



1 2 

log[r/(h- 1 Mpc)] 




log[r/(h" 



Figure 2. Real-space correlation functions of clusters extracted from the SCDM07 output and 03CDM at z=0 using either (F)OF 
or (H)OP algorithms. The two left panels show the correlation functions for clusters with masses greater than the specified threshold. 

— 1/3 

The numbers in the parenthesis are the D c = n c 1 values for the cluster samples. The two right panels show correlation functions for 
samples with the cluster number density given the specified value of D c . 



applying a fixed lower mass cutoff. We also generated clus- 
ter catalogs based on the ordered list with a specific number 
density of clusters (labelled by the corresponding value of 
AO- 

For each cluster catalog, we compute the real space two- 
point correlation function using the direct estimator: 



Mr) 



N P (r) 
n 2 c V(SV) 



(2) 



where N p (r) is the number of cluster pairs in the radial bin 
of volume SV centered at r, n c is the mean space density of 
the cluster catalog and V is the volume of the simulation. 
We use all the clusters in our catalogs, taking advantage of 
the periodic boundary conditions. 

The la error bars for the correlations are estimated 
using the formula 



3 1 
2 



■(i + £«(»■)), 



(3) 



where N cc is the number of distinct cluster pairs in the radial 
bin at r. We have increased the size of the Poisson error 
bars by 50% because these errors do not take into account 
clustering and so are likely to underestimate the true errors 
(Croft & Efstathiou 1994; Croft et al. 1997). 

The correlation functions are not well described by a 
single power-law over the entire range of pair separations 
sampled. To estimate the correlation length, we fit a func- 
tional form 



(4) 



over the range 4.5/i _1 Mpc < r < 25/i -1 Mpc, which brack- 
ets the point where £ cc = 1. We estimate the value of Ro by 
both fixing the value of 7 in the above equation to —1.8 (see 



© 1998 RAS, MNRAS 000, 



6 F.Governato et al. 




PS , 



1Q13 1Q14 | fl 15 1Q U 

M/h-'Mg 



Figure 1. Differential mass function (number density per unit 
mass) of groups and clusters extracted from the z = 03CDM 
simulation volume using FOF (continous line) and HOP (dot- 
dashed line) algorithms. The dotted line shows the analytic Press- 
Schechter prediction for the mass function. 

equation [l]) and by allowing 7 to be a free parameter. Since, 
the fit is done over a restricted range in r, both schemes 
yield similar values of Ro- 

In Figure ^, we show real-space correlation func- 
tions of clusters in catalogs defined by two different lower 
mass thresholds (M cut = 4.6 x 10 13 /i _1 M Q and 2.5 x 
10 14 /i _1 Mq) and two different cluster abundance require- 
ments (D c = 15k" 1 Mpc and 40/i _1 Mpc). The clusters are 
extracted from the simulations using either the FOF or HOP 
algorithms. Figure 2 shows the results for clusters extracted 
from the SCDM07 output and the 03CDM simulation at 
z=0; the clustering trends of 04CDM and SCDMfO clusters 
are the same. 

At both low and high mass thresholds, the correlation 
functions of FOF and HOP clusters are virtually identical, 
especially in the range 4.5/1" 1 Mpc < r < 25/i _1 Mpc. The 
abundances of FOF and HOP clusters (or equivalently, their 
D c value) are also the same, as expected from results shown 
in Figure f . As the mass threshold increases, or the number 
density is decreased, the clustering amplitude increases, but 
the shape of the correlation function remains the same. This 
reflects the fact that massive, rare peaks tend to be more 
strongly clustered in all CDM models. This result is shown 
in Figure |[ 

Given the good match between the halo catalogs we will 
mainly discuss results for the FOF clusters. Unless specified, 
results for FOF clusters hold for the HOP clusters as well. 



3.1 Cosmology and Normalization of the Mass 
Power Spectrum 

The real-space z = correlation functions of FOF clus- 
ter samples from the various simulations are compared in 
Figure §. Considering the SCDMfO and SCDM07 as two 
z—0 outputs results in the model with the higher amplitude 
(SCDMfO) developing structure on group and cluster scales 
at an earlier epoch, having a higher density of very massive 
halos and a more strongly clustered mass density field at the 
present epoch. 

In spite of the above mentioned differences, the cluster 
correlations for the SCDM models with different normali- 
sations are virtually identical in shape and amplitude for 
cluster samples with both high as well as low mass thresh- 
olds. In the case of massive clusters, this has been previously 
noted by both Croft & Efstathiou (1994) and Eke et al. 
(1996). Since for the SCDM models, studying the changes 
(or lack thereof) in the correlation functions due to varia- 
tions in the normalization of the amplitude of the primordial 
density fluctuations is equivalent to studying the evolution 
of the clustering property as a function of time, we defer the 
discussion of the above-mentioned until §3.3. 

The correlations for the two open models are also very 
similar to each other, both in shape and amplitude. These 
two models differ not only in their values of fio an d h but 
also in the normalization of the amplitudes of the primordial 
mass fluctuations as defined by as- The two OCDM models 
do, however, have the same value of floh 2 . Since it is this 
parameter that defines the position of the peak in the CDM 
power spectrum characterizing the initial Gaussian random 
fluctuations in density field, it is perhaps not surprising that 
the cluster correlations are similar. 

In comparison to the cluster correlations in a critical 
universe, the OCDM cluster correlation functions have a sig- 
nificantly higher amplitude. This occurs because the peak in 
the power spectrum for the OCDM models is displaced to- 
wards larger scales and therefore, for a similar value of as, 
the OCDM models have more power on large scales than the 
SCDM models. 



3.2 Redshift Evolution 

In Figure [B|, we plot the present-day and z = 0.43, z — 
0.58 cluster correlations, in comoving coordinates, for two 
cluster samples defined as (M > 4.6 x 10 13 /i _1 Mq and 
M > 1.5 x 10 14 /i _1 M Q ) drawn from the SCDMfO and 
03CDM models respectively. 

Before we discuss the results, let us consider what is 
expected. Given a sample of halos with masses greater than 
some threshold M cu t, the correlation function of the halos 
can be related to that of the total mass distribution via the 
bias parameter: 

£cc(r;M > M cut ) = b 2 eff (M cut )£ pp (r). At a given epoch 
the bias parameter becomes larger as M cu t is raised, as we 
have already shown. For a fixed M cu t and a critical universe, 
the bias parameter is expected to decrease asymptotically 
to unity as a function of time (Tegmark & Peebles 1998) as 
the underlying mass distribution becomes more clustered. 
The time evolution of £cc(r;M > M cu t) depends on the 
competition between these two trends. 

Turning to Figure B, we note that for both low and high 
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Figure 3. Real-space correlation functions of clusters samples defined by either imposing different mass thresholds (top two panels) 
or by demanding that the sample clusters have some predefined number density (bottom two panels). The clusters have been extracted 
from the SCDM07 output and 03CDM at z=0 using the FOF algorithm. 



mass thresholds, there are no significant differences between 
the comoving correlation functions at z = (crg=l) and 
z = 0.43 for the SCDM clusters. This implies that over the 
mass and redshift ranges considered here, the rate of increase 
in the clustering of the mass distribution is closely matched 
by the rate at which the bias parameter decreases. 

For the 03CDM model, the correlation function at 
the earlier epoch has a slightly higher amplitude for both 
low and high threshold samples. The comoving correlation 
length at z = 0.58 is a factor of 1.1—1.2 greater. 

To summarize, the comoving group / cluster correlation 
functions are either constant or change very little over the 
redshift range < z < 0.5 and in proper coordinates, the 
group /cluster correlation length decreases with increasing 
redshift over the redshift range studied. In the SCDM case, 
this decrease is given by Ro oc (1 + z) -1 and in the 03CDM 
model, by R oc (1 + z) -0 ' 86 . 



Table 1. The bias parameter for SCDM07 cluster samples. £>i 
is computed using the standard PS mass function. The second 
column gives the mass cut in units of the characteristic mass, 
M* = 4 X 10 13 h _1 Mq. The final column gives the Lagrangian 
radius of the halo, which is the smallest separation where the 
assumtpions in the calculations are valid. 



Mcuth' 1 Mq 


Mcut/M* 


6i 


r L h- x Mpc 


5.7c+14 


16.2 


3.3 


7.9 


2.8c+14 


7.0 


2.4 


(5.2 


2.1c+14 


5.2 


2.2 


5.7 


7.0c+13 


1.7 


1.5 


3.9 



© 1998 RAS, MNRAS 000, 



8 F.Governato et al. 




M > 1.5xl0 14 h- 1 M Q - 
• SCDM07 (FOF) I 




1 2 

log[r/(h _1 Mpc)] 

Figure 4. Real-space z = cluster correlation functions 
extracted from our simulations (SCDM07, SCDM10, 03CDM, 
04CDM) using the FOF algorithm. 



Table 2. The bias parameter for 03CDM cluster samples, as in 
Table [i] In this case M, = 1.4 X 10 13 /i _1 M 



Mcnth' 1 Mq 


Mcut/M* 


&l 


ri.hr 1 Mpc 


2.7c+14 


19.3 


2.7 


9.2 


1.9c+14 


13.6 


2.4 


8.2 


l.lc+14 


7.8 


2.0 


6.8 


9.0c+13 


6.4 


1.9 


6.1 


1.8c+13 


1.3 


1.2 


3.7 
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Figure 5. Correlation functions of SCDM10 and 03CDM clus- 
ters computed at two different epochs and plotted in comoving 
coordinates. 



3.3 Comparison with analytic calculations 

To date, most studies of cluster correlations have utilized 
numerical simulations. Such numerical simulations are very- 
expensive to generate, a constraint that renders a systematic 
exploration of different cosmological models impractical; it 
also makes it rather difficult to explore and identify the gen- 
eral physical mechanisms underlying the clustering proper- 
ties of clusters and group halos viz a viz that of the mass dis- 
tribution. Consequently, various authors (e.g. Kaiser 1984; 
Bardeen et al. 1986; Mann et al. 1993; Mo & White 1996; 
Mo, Jing & White 1996, Catelan et al. 1998) have devel- 
oped analytic schemes to compute the cluster correlation 
function. 

The first method to compute cluster correlation func- 
tions analytically that we discuss is based on the Press- 
Schecter formalism and its extensions. This was originally 
developed by Cole & Kaiser (1989) and Mo & White (1996) 
to derive a model for the spatial correlation of dark mat- 
ter halos in hierarchical models. The calculation consists of 
three steps (see Baugh et al. 1998). 

• Compute the nonlinear power spectrum for the cos- 
mology and as in question using the transformation of the 
linear power spectrum suggested by Peacock and Dodds 
(1996). 

• Calculate an effective bias parameter, fe e // for the 
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dark matter halos above the specified mass cut as outlined 
by Mo & White (1996). 

• Fourier transform the nonlinear power spectrum to get 
the nonlinear correlation function of the mass, then multiply 
by the square of the halo bias factor, to get the real-space, 
nonlinear, cluster correlation function: £cc{r) = ^e//Cpp( r )- 

The cluster correlation function thus computed has 
been tested, against N-body results by Mo & White (1996) 
and Mo, Jing & White (1996) and is found to hold even 
in the mildly nonlinear regime where f(r) > 1 as long as 
r > Tl where tl = (SM/4-Kpo) 1 ^ 3 is the Lagrangian radius 
of the dark matter halos (rz, ~ 10/i -1 Mpc for rich clusters 
of galaxies) and po is the present mean density. Recently 
Jing (1998) has shown that the Mo & White formula sys- 
tematically underpredicts the bias of low mass halos, but it 
is in good agreement with numerical simulations in the mass 
range considered here. 

The bias parameter for a dark matter halo that con- 
tains a single galaxy is given by the formula derived by Mo 
& White (1996) and was written down for any redshift in 
Baugh et al. (1998): 



b(M,z) = l+ T 
8 C 



a(M)D{ 



- 1 



(•») 



Here D(z) is the linear growth factor, normalized to 
D(z = 0) = 1, a(M) is the rms linear density fluctuation 
at z — and S c is the extrapolated linear overdensity for 
collapse at redshift z. This gives the bias factor for the halo 
when the clustering is measured at the same epoch that the 
halo is identified. 

For a sample of halos with different masses, the effective 
bias is given by 



Oeff 



(*) = 



/ N(M, z) b{M, z) dM 
J N(M, z) dM 



(6) 



where N(M, z) dM is the number density of halos with 
mass M in the sample. For the cluster samples that we 
have constructed, N(M, z) can either be set equal to the 
Press-Schechter mass function (with 5 C set to the canonical 
value defining collapse for the cosmology under considera- 
tion) or to the cluster mass function computed directly from 
the cluster catalogs. However the results shown here are al- 
most insensitive to this choice. 

The effective bias parameters for samples whose corre- 
lation functions are plotted in Figure |^ are given in Tables 
[l] (SCDM07) and | (03CDM). The mass cuts applied corre- 
spond to halos of different rarity in the two cosmologies; this 
is quantified by comparing the mass cut to the characteris- 
tic mass M, , which is defined later in §4. All the mass cuts 
considered correspond to objects that are greater than M*, 
and so these halos are biased tracers of the dark matter dis- 
tribution (Mo & White 1996). The cluster sample with the 
highest mass cut for both 03CDM and SCDM07 is predicted 
to have a correlation function that is ~ 10 times higher than 
that of the dark matter. 

The PS-based analytic correlation functions are shown 
in Figure ^ as solid curves. There is little difference between 
the correlation functions computed using the standard PS 
mass function and the numerical mass function discussed 
in §4. The analytic correlations are in excellent agreement 
with our numerical correlation functions. The agreement be- 



tween the numerical and analytic results is further confirmed 
by the match between the analytic and numerical Ro~D c 
curves. The analytic Ro~D c curve is plotted in Figures [?] 
and |8] as the light solid curve. 

We consider next the scheme developed by Mann et al. 
(1993). This is based on the method devised by Couchman 
& Bond (1988; 1989) that combines the theory of the statis- 
tics of peaks in Gaussian random fields with the dynamical 
evolution of the cosmological density field. 

In this scheme, the time evolution of the density field 
is followed using the Zel'dovich approximation (Zeldovich 
1970). At the epoch of interest, a particular class of objects 
is defined by the pair R s and 5 C . These are, respectively, the 
smoothing scale that is applied to the cosmological density 
field and the linearly extrapolated amplitude of the density 
fluctuations at the time of collapse. Mann et al. (1993) set 
the values of these two parameters by choosing an appro- 
priate value for S c (S c = 1.686, corresponds to collapse of 
spherical density perturbations in an Q = 1 universe) and 
then adjusting R s until the number density of peaks with 
overdensities greater than 8 C corresponds to number density 
of objects under consideration. Full details can be found in 
Mann et al. (1993). 

In Figure [], we plot the correlation functions of some 
of our cluster samples and show the corresponding analytic 
peaks-based correlation functions computed assuming S c — 
1 and S c = 1-7 (dashed curves) according to Mann et al. 
(1993)'s prescription. 

As noted by Mann et al. (1993), the correlation func- 
tions computed using 8 C = 1.7 consistently overestimate the 
correlation amplitudes on all scales of interest. The correla- 
tion functions for 8 C — 1 are in excellent agreement with the 
numerical results for D c < 30/i _1 Mpc (however, a value of 
8 C ~ 1 is rather unphysical). For larger values of D c , the an- 
alytic results tend to overestimate the correlations, with the 
discrepancy first becoming obvious on small scales and then 
propagating out to larger scales as D c continues to increase. 
For a given D c , the discrepancy is more severe for SCDM07 
clusters than for 03CDM clusters. 

The tendency for the peak scheme to overestimate the 
correlations on scales where £,cc(r) > 1 by a margin that 
grows larger with increasing D c suggests that the peaks- 
based method will overestimate the correlation lengths of 
samples with large values of D c . Both Mann et al. (1993) 
and Watanabe et al. (1994) have computed the Ro~D c rela- 
tion predicted by the peaks method. Comparing their Ro~D c 
curves for 03CDM-like models with our numerical results 
we find that for cluster samples with D c ~ 80ft" 1 Mpc, 
the correlation length predicted by the peaks-based analytic 
scheme is Ro ~ 28/i _1 Mpc whereas the simulation result 
is Ro ~ 22/i _1 Mpc, very close to the PS based prediction. 
The breakdown in the peak scheme may be due to several 
factors. The simplest possibility is that that the method uses 
the number density of clusters in a sample as a constraint 
rather than some physical properties of the clusters such as 
their masses. Other possible causes of the breakdown are: 
the manner in which the different filtering scales are chosen, 
the simplistic nature of the prescription defining the relation 
between peaks in the smoothed density field and "clusters" , 
and the requirement that the large-separation asymptotic 
limit of the statistical contribution to the cluster correla- 
tion function matches the statistical peak-to-peak corre- 
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Figure 6. Analytic correlation functions compared against our numerical results. The two dashed curves are the peaks-based correlation 
functions computed according to the prescription of Mann et al. (1993). The curve with the higher correlation amplitudes on small 
scales corresponds to <5 C = 1.7. The other curve corresponds to <5 C = 1.. The two solid curves are the PS-based correlation functions 
computed as described in §3.5. The correlation functions are computed assuming either the standard PS mass function or the numerical 
cluster mass function (see §4). The two are very similar. 



lation function (Mann 1998, private communications). The 
latter two tend to magnify any minor discrepancy caused by 
any of the other factors. 

3.4 Correlation Length and the Cluster 
Abundance 

In Figures [7] and ^ we plot correlation length (Ro) as a func- 
tion of cluster abundance in terms of D c , the mean cluster 
separation, for the FOF clusters in the SCDM07 output and 
03CDM model at z=0. For comparison, we also show the 
scaling relation (Equation jij); the numerical results of Bah- 
call & Cen (1992) and Croft & Efstathiou (1994) the obser- 
vational data for R> 0, R> 1, R> 2 Abell clusters (open 
triangles) from Bahcall & Soneira (1983) and Peacock & 
West (1992) and the data for the APM clusters (open cir- 



cles) given by Dalton et al. (1992) and Croft et al. (1997). 
In neither of the SCM07 nor 03CDM models is the R - 
D c relation for clusters consistent with the scaling relation 
Ro = 0AD c . 

However, it should be noted that the numerical results 
show the Ro-D c relation for the real-space correlation func- 
tion whereas the correlation lengths for the observed clus- 
ters are derived from redshift-space correlation functions. 
Redshift-space correlation lenghts are generally larger than 
their real-space correlation length counterparts. For exam- 
ple, in the case of their low-density spatially flat CDM 
model, Croft & Efstathiou found that their real-space Ro-D c 
relation saturates at Ro Ri 15/i _1 Mpc for large values of D c , 
whereas the redshift-space Rq-D c saturates at Ro ~ 21/i _1 
Mpc. An increase of this kind, however, is not sufficient to 
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Figure 7. Cluster correlation length as a function of D c , the mean cluster separation, for clusters extracted from the SCDM07 output 
using FOF (heavy solid curve). The error bars show the la errors in Ro derived from fitting the correlation functions with a —1.8 
power-law as described in the text. The solid line corresponds to the scaling relation between Rq and D c (equation |l|) advocated by 
Bahcall & West (1992). The dot-dashed line shows the Rq-D c that Bahcall & Cen (1992) derived from their numerical study. The 
short-dashed curve corresponds to the numerical results of Croft & Efstathiou (1994). In addition, the open triangles show the results 
for R> 0, R> 1 and R> 2 Abell clusters (Bahcall & Soneira 1983; Peacock & West 1992) and the open circles show the results for APM 
clusters (Dalton et al. 1992; Croft et al. 1997). The light solid curve is the Rq-D c relation derived from analytic PS-based correlation 
functions computed according to the prescription in Baugh et al. (1998) — see §3.5. 



bring our numerical Rq-D c relation into agreement with the 
scaling relation of equation [j]. 

Consistently with all previous findings, the Ro-D c curve 
for clusters extracted from the SCDM universe (Figure 
does not match either the APM or Abell results. On the 
other hand, the results of our 03CDM model are in good 
agreement with the APM and richness R> 0, R> 1 Abell 
cluster data, even if the effect of redshift distorsions increas- 
ing the length scale Ro a few Mpc were included. 

The seriously discrepant datapoint is for R> 2 Abell 
clusters. If this measurement is correct, it suggests that 
clustering on very large scales may have been modulated 
by non-Gaussian processes (see Mann et al. 1993; Croft & 
Efstathiou 1994) as it is very difficult to conceive of a Gaus- 
sian model that can produce the requisite clustering at these 
scales. It would also imply that the very rich APM clus- 
ters with their comparitively large correlation lengths are 
not really rich or massive systems but rather are systems 
comparable to R> 1 Abell clusters whose number densities 
have been biased downward by the cluster identification al- 
gorithm. The agreement between our analytic results and 



our numerical results for clusters in the 03CDM model and 
the APM results leads us to believe that it is the R> 2 Abell 
result that is most likely incorrect, biased upward by the in- 
homogeneities and contamination due to projection effects 
in the Abell catalog as argued by Sutherland (1988) and 
Sutherland & Efstathiou (1991). 

In comparing our numerical results for SCDM07 to 
those of Bahcall & Cen (1992), we find that for 25/i -1 Mpc 
< D c < 40/1" 1 Mpc, our Ro-D c results are consistent with 
theirs. For D c > 40/i -1 Mpc, our curve rises less steeply 
than that of Bahcall & Cen and appears to saturate for 
D c > 50h~ 1 Mpc. From analytic results (light solid curve), 
which we discuss further in the next subsection, we expect 
the Ro-D c curve to continue to rise but much more gently 
than the Bahcall-Cen result. Since both we and Bahcall & 
Cen (1992) used the FOF algorithm to identify clusters in 
the simulations, the cluster selection algorithm cannot be 
responsible for the differences. Furthermore, our correlation 
lengths were determined in the same way as Bahcall & Cen 
(1992). 

Comparing our 03CDM results to those of Bahcall & 
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Figure 8. Same as Fig. |^ but for 03CDM model. 



Cen's (1992) low-f2 models, we find that the two are in good 
agreement for D c < 35ft _1 Mpc and also in a rough agree- 
ment with the scaling relation. However as the cluster abun- 
dance decreases and D c increases, the correlation lengths of 
our cluster samples do not increase as quickly. Our numeri- 
cal results at large values of D c are in good agreement with 
those derived analytically to the scales probed by our simu- 
lations, (see §3.3) 

In comparing our as = 0.7 result to Croft & Efstathiou's 
(1994) r c = 1.5ft," 1 Mpc a 8 = 0.59 SCDM model, we once 
again find smaller correlation lengths for D c < 30ft" 1 Mpc. 
For higher values of D c , the Croft & Efstathiou (1994)'s re- 
sults are consistent with ours in spite of the fact that wc 
have used FOF to identify the clusters and Croft & Efs- 
tathiou (1994) results are based on a very different scheme. 
Comparing our 03CDM results to those of Croft & Efs- 
tathiou's (1994) 1.5ft" 1 Mpc, a 8 = 1.0 low-Q spatially fiat 
CDM model, we find that within the uncertainties in the two 
curves, they are in excellent agreement with each other. The 
flattening in Croft & Efstathiou's curve for D c > 50/i _1 Mpc 
(at Ro ~ 15ft -1 Mpc) is not real. As indicated by both our 
numerical and analytic results, the correlation continues to 
rise, albeit gently, reaching Ro ~ 22ft -1 Mpc at D c — 80ft _1 
Mpc and is still rising. The flattening trend is likely an 
artifact of the finite simulation volume or even the poor 
mass/force resolution. 



4 THE CLUSTER MASS FUNCTION 

According to the analytic PS formalism, the comoving num- 
ber density of dark matter halos of mass M in the interval 
dM is 



N(M) = J 2 --P ^ 



TV M 2 



diner 



dlnM 



exp 



SiD~ 



2a 2 



(7) 



where p is the comoving density of the Universe and a(M) 
is the linearly extrapolated present-day rms density fluctu- 
ation in spheres containing a mean mass M. The redshift 
evolution of N(M) is controlled by the density threshold 
for collapse, S c /D(z), where D(z) is the linear growth fac- 
tor normalized to unity at z — (Peebles 1993) and S c is 
the linearly evolved density contrast of fluctuations that are 
virializing at z — 0. The growth factor, D(z), depends on Qq 
and A whereas S c has only a weak dependence on Qq. For 
spherical density fluctuations, S c — 1.686 for fi = 1 and 
1.65 for fi = 0.3. 

The PS description of structure formation in the Uni- 
verse leads naturally to the definition of a characteristic 
mass M»(«) such that 



a{M*)D{z) = 5 C . 



(8) 



M* (z) is then the characteristic mass of halos that are virial- 
izing at redshift z. Its evolution tracks the manner in which 
structure forms. In bottom-up hierarchical clustering mod- 
els, such as CDM models, M,(z) increases as a function of 
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time as lower mass structures are incorporated into progres- 
sively more massive halos. In a critical universe, the growth 
factor evolves as D(z) = (1 + z)^ 1 and to first order, this 
implies a strong evolution in M*. In an open or a flat, low-f2o 
universe, D(z) ceases to evolve as strongly, and the evolu- 
tion of the characteristic mass is greatly suppressed, once 
fi(z) deviates significantly from unity. Hence, the evolution 
of the dark halo mass distribution is also greatly suppressed. 
A clear detection of the presence or absence of strong dy- 
namical evolution in the cluster population can be used to 
put stringent limits on the underlying cosmology. 

The actual value and the details of the evolution of 
M*, and therefore of the mass distribution especially for 
M > M», depends sensitively on 8 C . The standard practice 
is to use the value of 5 C for collapse of spherical perturba- 
tions. Typical perturbations in CDM models, however, are 
not spherical and therefore, the actual value of S c will dif- 
fer from the spherical value. Moreover, Heavens & Peacock 
(1986), argue that S c is likely to be lower (1 < 5 a < 1.68) 
because typical proto-structures in a Gaussian random field 
tend to be triaxial. A lower (higher) value of 8 C results in 
more (fewer) high mass objects. A detailed discussion of S c 
and the asphericity of the density perturbations is given by 
Monaco (1995, 1998), who finds that when the assumption 
of spherical collapse is relaxed, 8 C becomes a function of the 
local shape of the perturbation spectrum. In most cosmo- 
logical models, the power spectrum of the primordial per- 
turbations over the scales of interest deviates, albeit gently, 
from a simple power law shape and it becomes debatable 
whether a constant value of 8 C is a fair description of the 
evolving cosmic mass function at all masses and redshifts. 
Here we allow the collapse threshold to be a free parameter 
depending on redshift, calibrating at each epoch using the 
high mass end of the mass distribution. As discussed in §1 
and §2.3, the halo mass function derived using the PS for- 
malism is a measure of the abundance of collapsed, distinct, 
halos characterized by their virial radius and mass. Conse- 
quently, it is appropriate to use catalogs generated using the 
FOF and HOP cluster finding algorithms. 

4.1 Computing the halo mass function 

We construct the differential mass function by sorting the 
halos according to their masses in bins of size Alog(M) = 
0.1. We have verified that our results are insensitive to this 
choice of bin size. Due to the large size of our simulations, we 
are able to study also the differential cluster mass function 
instead of the cumulative distribution, as is usually done. 
This means that the individual bins are independent and 
the results more robust. We estimate the uncertainty in the 
number of objects in each bin using Poisson statistics. It 
is useful to remember that the SCDM run can be rescaled 
freely to a different erg normalization, The redshift of each 
given output is then rescaled to a redshift z': 1 + z' = (1 + 
z)/a s . 

4.2 Effects of numerical resolution and cosmic 
variance 

To study the effects of degrading the numerical resolution, 
we examined a lower resolution run of our SCDM07 volume. 



This run used the same phases, 3 million particles, a soften- 
ing length of 160/i -1 kpc and a third of the timesteps used in 
our fiducial run. Both force and spatial resolution are there- 
fore significantly poorer. Halos represented by 128 particles 
in our fiducial run have only 8 particles in the low resolu- 
tion simulation. While there is good agreement at the high 
mass end, the low mass end of the mass function is severely 
affected by the poorer resolution, showing a significant de- 
crease in the number density of halos below 3 x 10 14 /i _1 Mq. 
This test indicates that at least 30 particles are needed to 
correctly assign a mass to individual halos. Our choice to 
include in our analyses only halos with N> 64 is then a 
conservative one. 

Finally, we also explored the effects of cosmic vari- 
ance on the cluster number density. We divided our z = 
SCDM07 volume into several subvolumes and measured the 
local 8 C over the same mass range as we did for the whole 
volume. As expected, cosmic variance produces a scatter in 
S c when measured in smaller volumes. However the scatter 
is not significant for volumes as little as 1/8-th of the origi- 
nal simulation volume (i.e. cubes with 250h~ 1 Mpc per side): 
we find S c = 1.68 ± 0.02. This is close to the error associ- 
ated with a single measurement and similar to the value we 
get from the whole volume, suggesting that the value of S c 
for a given halo finder has (almost) converged when the full 
volume of the simulation is considered. 

4.3 Press-Schechter Predictions vs Numerical 
Results 

In Figures ^| and [l^, we show the differential mass func- 
tions from the SCDM and 03CDM simulations, respectively. 
The corresponding PS curves, computed using the canonical 
value of S c for spherical perturbations, are also shown. We 
only show the FOF results. Generally, the HOP and FOF 
results are very similar, with HOP having a slightly larger 
number of massive clusters. 

In the case of the SCDM model, we find that the shape 
of the differential PS mass function is roughly consistent 
with the shape of the numerical mass function only at erg = 1 
At lower as (or alternatively, higher z) the PS mass func- 
tion underestimates the number density of rich clusters in 
the simulation. The excess at the high mass end is about a 
factor of a few in number density per mass bin. At as 0.7 or 
larger the PS approach overestimates the number of small 
halos (M < lO 14 ^ 1 M Q ). The deficit of low mass halos 
(which we will only touch upon briefly here) has been well- 
documented in numerical works by Carlberg & Couchman 
(1988) and more recently, by Lacey & Cole (1994), Gross et 
al. (1998) and Somerville et al. (1998). This deficit arises 
independently of the choice of algorithm used to define the 
halos in the simulations (see Figure [l]) and has been as- 
sociated with merger events not accounted for by the PS 
formalism (see Cavaliere & Menci 1997 and Monaco 1997). 
The fact that the two halo finders agree extremely well in 
this regime makes the result very robust. 

Apart from the above-mentioned deficit, most other 
studies that have tested the analytic PS mass function 
against numerical results have reported a good agreement 
between the two (e.g. Eke et al. 1996) with the notable ex- 
ceptions of Bertschinger & Jain (1994) and Sommerville et 
al. (1998), who found a systematic excess of massive halos 
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Figure 9. The N-body mass distribution for the SCDM run — shown as points — at four different outputs as well as the standard PS 
mass function (computed using S c = 1.686) — shown as lines — for the same four outputs: erg = 1 (filled circles, solid line), trg = 0.7 
(circles, dotted line), its = 0.47 (squares, dashed line), trg = 0.35 ( triangles, long dashed line). The error bars correspond to la Poisson 
errors. 
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Figure 10. The same as Figure El but for 03CDM model at z=0,0.58,l. 
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in numerical simulations as compared to the PS prediction. 
These numerical results, however, are derived from simula- 
tions that have lower resolution and probe smaller cosmolog- 
ical volumes than our simulations. These simulations, con- 
sequently, contain only a small number of massive clusters 
and this, in conjunction with cosmic variance, has resulted 
in numerical mass functions with large uncertainties at the 
high mass end. Eke el al. (1996), found a good agreement 
with the PS mass function on the scale of 5 x 1O 14 /i _1 M0 
but with a large error bar. 

The large volumes used in our work allow us to confirm 
the good agreement between the numerical and the analyt- 
ical mass function found for high values of erg. The excess 
of massive clusters at lower erg /higher redshifts suggest that 
the cluster mass function for the SCDM model evolves more 
slowly than the PS mass function. At ag=0.35, the number 
density of Coma-like clusters exceeds the standard PS pre- 
diction by almost an order of magnitude. Allowing S c to be 
a free parameter at each epoch, we perform a % 2 fit of the 
PS formula to the numerical results assuming Poisson errors 
for each bin. Since observations tend to be biased towards 
the high-mass end, we include in our fit only clusters with 
temperature T>3keV, (the M-T relation is defined in §4.4) 
to allow the use of this formula in theoretical predictions for 
X-ray clusters. This formula will underestimate slighlty the 
number of very (T>7keV) hot clusters for high erg models. 

For the SCDM cluster mass function at og=0.7 , the 
best-fit S c coincides with the canonical value of 1.686 ! For 
z > and for our FOF selected halos, the best-fit value of 
S c as a function of redshift and erg is : 




Figure 11. The points show best fit values of S c (z) required 
to match the numerical mass function for clusters with T> 3 
keV at different redshifts. The squares show the results for the 
SCDM07 model and the triangles show the results for 03CDM 
model. The horizontal light solid curve and the nearly horizontal 
dashed lines show the canonical values of 8 c (z) for the SCDM07 
and the 03CDM models respectively. The lines across the points 
are the power-law interpolation to the points (see text). Bars are 
3 a errors. 



8c(z) = 1.686 



07 

0"8 



(1 + *) 



(9) 



as shown in Figure [ll], whereas the canonical value of 
S c is a constant, (in Figure pd| errorbars are 3 a errors). For 
HOP selected halos S c is offset toward even lower values, 
i.e. toward a larger cluster excess: 8 C — 1.6 at z=0. (for the 
SCDM model at ag = 0.7). However, a similar evolution of 
S c (z) is found, with a (1 + z) -0 ' 1 dependence. In principle 
the fitting procedure should keep into account that the mass 
associated with a given temperature depends on z and that 
the same output can be associated with different z depending 
on erg. However the fitting formula is not affected by this for 
values of present day erg between 1 and 0.5 and eq. ^ can be 
used safely. 

In the 03CDM model case , the PS mass function is in 
fair agreement with the numerical mass function, especially 
at low redshifts. At z = 0, the high mass end of the numer- 
ical halo mass function agrees within 2a with the analytical 
curve computed using the canonical (spherical) value of S c 
= 1.651. The uncertainties in the number densities of mas- 
sive clusters are slightly larger in the 03CDM case because 
of the smaller physical volume/h igher Ho of the simulation. 
As shown in Figure O, the best- fit S c (z) for FOF halos can 
be well-approximated as 



5c{z) = 1.775(1 + z) 



(10) 



HOP results shows an even weaker evolution in z. As the 
result is more significant in the SCDM case we will mainly 
focus on the analysis of results for the critical case. 

From these results however, it is not clear if the devia- 
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Figure 12. The points show the best-fit values 5 c (z) for SCDM07 
(continuous) and 03CDM (dashed) at z=0 (black dots) and z=l 
(i.e <rg = 0.35 for SCDM) (open squares) as the mass range over 
which the PS functional form is fitted to the numerical FOF mass 
function is varied. The abscissa corresponds to the lower mass 
threshold of the mass range over which the fit is demanded. 
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tion from the standard PS mass function is due to just one 
or rather both of the following effects: 

• at lower as and for a fixed temperature T we study 
more extreme clusters i.e. we look at a different region of 
the mass function, which maybe still be self similar, albeit 
different from the canonical PS. 

• the shape of the mass function evolves with time 
and/or depends on the power spectrum 

As discussed above, the best-fit S c (z) were determined 
by fitting the functional form of the PS mass function to the 
numerical results for clusters with T> 3 keV, i.e. on mass 
scales M > M*, where M» corresponds to 4 x 10 13 h~ 1 M Q 
for SCDM07 and 1.5 x 10 13 h~ x M Q for 03CDM. It is of 
interest to relax this constraint and explore how the S c varies 
as the minimum mass of the halos included in the fit is 
lowered and approaches the mass of the smallest halos (64 
particles) in our catalog. We carried out the above exercise 
at two different epochs/normalization and the results are 
shown in Figure hi The value of 5 C changes dramatically as 
the mass range over which the fit is carried out moves toward 
smaller masses (the fit is dominated by the smaller mass 
bins as they contain most of the halos used in the fitting). 
The trend for both SCDM and 03CDM is for 8 C to become 
larger as the mass threshold is lowered. This is precisely 
what one expects given the deficit of low-mass halos in the 
simulations. This shows that the shape of the N-body mass 
function differs from the PS prediction, and the exact value 
of 5 C is a function of the mass interval considered. 

This plot shows also that, at a given mass, 8 C is a func- 
tion of redshift. To prove that the shape of the mass func- 
tion evolves with time we take advantage of the fact that 
within the PS framework the cumulative fraction of mass in 
collapsed halos is invariant when plotted vs the variance of 
the density field at a given mass/lenght scale. I.e. at a given 
value of a (which corresponds to different mass/lenght scales 
depending on cosmology and z) the fraction of mass in col- 

Ssed objects is always the same. This is shown in Figure 

It is interesting to interpret this change in terms of dif- 
ferent power spectra. SCDM models can be rescaled to mod- 
els with different V' (or rCDM models) rescaling the box by 
r/r' and choosing as the final output the one with the cor- 
rect present-day normalization at the scale corresponding 
to the 8h _1 Mpc scale. The sequence of outputs with de- 
creasing as and increasing excess of massive objects can be 
reinterpreted as a sequence of models with smaller F param- 
eter and the same normalization, revealing the dependence 
of the shape of the mass function from the power spectrum. 
An excess of massive halos for a more negative local spec- 
tral index (as the case with larger F) had been predicted by 
Monaco (1995). 

This rescaling allowed to compare our mass functions 
with that obtained from the so called "Hubble Volume Sim- 
ulation" (HVS) (Colberg et al. 1998), a rCDM with r=0.21 
and normalization of as = 0.6 Their final output can be 
rescaled to a SCDM model with a a 8 = 0.285. The HVS 
mass function lies nicely along our sequence of mass func- 
tions (Cole & Jenkins 1998, private communication) showing 
a slightly larger excess of halos compared to our as = 0.35 
output (our output with the lowest as)- 

The two simulations used completely independent soft- 
ware to generate the initial conditions, evolve the density 



field and analyze the data. The agreement found is quite 
satysfying and shows that the results of both simulations 
are free of hidden systematic effects. 

We can conclude that for critical CDM models the nu- 
merical mass function compared to the PS analytical predic- 
tion has an excess of halos for M^> M» and a deficiency 
for masses M -C M*. The excess at large masses is larger 
for models with smaller F (or more negative local spectral 
index). 

These deviations for the canonical predictions are sig- 
nificant and cosmological tests based on the number density 
of a particular class of objects need to use the appropriate 
value of S c to make robust predictions. Our fitting formula 
eq. ^ can easily be modified for other critical models with a 
different shape parameter. 

4.4 Effects on the Cluster Temperature 
distribution 

X-ray observations allow one to directly determine the clus- 
ter temperature or the cluster luminosity function. Of the 
two, the temperature of the intracluster gas is thought to 
be the more robust measure of the depth of the cluster po- 
tential well and to a good approximation, is expected to 
be very strongly correlated with the cluster mass. In this 
section, we examine the impact on the cluster temperature 
function of the excess of massive clusters found in the simu- 
lation as compared with the predictions of the standard PS 
mass function. 

We follow Eke, Cole & Frenk (1996) and adopt the fol- 
lowing simple relation for estimating the temperature of the 
intracluster medium of a cluster of mass M : 

fc T = 7.75M^(l + ,)(^) 1/3 . (11) 

Ac is the average density contrast at virialization with re- 
spect to the critical background density. This relation as- 
sumes that the intracluster medium is isothermal. M15 is 
mass in units of 10 15 Mgh -1 

The cumulative temperature function iV(> T) consid- 
ering the SCDM07 model as the present time is shown in 
Figure [li]. The corresponding cumulative temperature func- 
tion based on the standard PS mass function is shown by 
the dashed lines. The bottom panel shows the ratio between 
these two mass functions. At z = 0, the two are very similar. 
At z ~ 1, however, the number density of clusters with tem- 
peratures kT > 7 keV (the temperature of rich, Coma-like 
clusters) is more than a factor of 3 greater than the canonical 
PS prediction and in the case of exceptionally hot clusters 
(kT > 10 keV), the discrepancy is an order of magnitude. 
This discrepancy increases for a universe with a lower, more 
realistic normalization. 

One implication of this is that the cumulative temper- 
ature function obtained from the simulation evolves more 
slowly over the redshift range < z < 1 than PS theory pre- 
dicts. Since we choose a fixed temperature range, at higher z 
we are observing more extreme (massive) clusters, i.e. a re- 
gion of the mass function with a stronger excess compared 
with the anlytical formula. In principle it may be more dif- 
ficult to challenge a high £1 model on the basis that the 
observed cluster temperature function does not vary signif- 
icantly over the redshift range < z < 0.3, especially if 
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Figure 13. The plot shows the cumulative fraction of mass in collapsed halos as a function of the variance of the density field. The PS 
prediction (invariant with redshift and cosmology) is the thick continuous line. The SCDM results are shown at diffferent outputs: erg = 1: 
dot-long dashed; erg = 0.7: dot-short dashed; erg = 0.47: dottted line; erg = 0.35: dashed. The numerical mass function is obviously not 
invariant and shows evidence of evolution. 
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Figure 14. Upper panel: PS (thick lines) vs numerical (thin lines) cumulative temperature function N(>T) for the SCDM07 model at 
2 = (continuous line), z = 0.5 (dotted), and z = 1.0 (dashed). Bottom panel: shows the ratio of clusters above a given temperature 
between the numerically-based temperature function and the PS-based temperature function at z = 0.5 (dotted), and z = 1.0 (dashed). 
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this effect is coupied with a significant scatter in the M-T 
reiation. The higher number density of hot ciusters in the 
simulations also suggests that attempts to estimate as by 
fitting the canonical PS-based temperature function may in 
principle lead to systematically high results. 

Following Eke et al. (1996, 1998) we reanalysed the 
low redshift cluster data of Henry & Arnaud (1991) and the 
new data in Henry et al. (1997) to produce a cumulative 
cluster temperature function. We then compared it with PS 
prediction modified using our inferred values for S c (z). This 
leads to a small revision of the amplitude of fluctuations 
inferred from the local cluster abundance for Qo = 1, to a 
value of as = 0.5 ± 0.04 if FOF halos are used. The use of 
HOP halos results suggests as = 0.48 ± 0.04. 

With a as= 0.5 normalization and including the rescal- 
ing of S c (eqVA) the excess of hot clusters with kT > 7 ke V 
compared to the standard PS prediction amounts almost to a 
factor of 10 at redshift of 1. Clearly, claims that rule out crit- 
ical density models on the basis of the detection of a single 
massive cluster at high z (Donahue et al. 1998) must then 
be taken with caution. However, the discrepancy between 
the temperature function predicted in a critical density uni- 
verse and that observed at z = 0.33 is reduced by a modest 
amount using the modified Press-Schechter scheme. The dis- 
crepancy there is still large enough to rule out Oo = 1, unless 
there are significant differences in the relation between mass 
and temperature for clusters at high and low redshift that 
make equation hll invalid. 



5 CONCLUSIONS AND DISCUSSION 

We have analyzed parallel N-body simulations of three 
CDM models: a critical density model (h=0.5) and two open 
models (Qo =0.3 and 0.4). These three models, span a large 
range of different properties (models COBE and cluster nor- 
malized, different amounts of large scale structure, high and 
low Ho) to cover the wide range of cosmological models 
presently considered. We simulated very large volumes - 
500h~ 1 Mpc per side - and used 47 million particles for each 
run. Having good mass, force and spatial resolution, these 
datasets allow a robust determination of two very important 
quantities for cosmological studies: the correlation function 
of galaxy clusters and the shape and evolution of the cosmic 
mass function on cluster mass scales. In particular, we have 
focused on how results from simulations compare with pre- 
dictions from analytical methods based on the PS approach. 

• Halo finders We have assessed if cluster correlations 
and mass functions are affected by the biases introduced by 
using different halo finders. This is an important step both 
for comparing the simulation results with theoretical predic- 
tions as well as with observational data. We used two halo 
finders available in the public domain: FOF and HOP. We 
found that the two halo finders, once set to select virialized 
structures based on a local overdensity criterion, produced 
very similar halo catalogs. While a small mass offset is de- 
tectable, (of the order of 5-7% for massive clusters) our con- 
clusions do not depend on the particular choice of the halo 
finder. In most of our work we conservatively used FOF, 
which gives results closer to the analytical approach. 

• Cluster correlation function We were able to deter- 
mine the clustering properties of halos spanning nearly two 



orders of magnitude in mass, ranging from groups and poor 
clusters to very rich massive clusters. Our analysis has shown 
that the comoving correlation functions, for a wide range of 
group/cluster masses, do not change significantly over the 
redshift range < z < 0.5. Of the two analytic schemes we 
used (peaks-based and PS-based) for computing cluster cor- 
relations functions, we found that the correlation functions 
derived using the PS-based scheme (Mo, Jing & White 1996, 
Baugh et al. 1998) were in excellent agreement with the nu- 
merical correlation functions. 

• Cluster Correlation Length and Number Density We 
compared our results for the cluster Rq-D c relation in dif- 
ferent models against both observations and results of pre- 
vious numerical studies. Firstly, we have confirmed that the 
SCDM model does not have sufficient large scale power to 
account for the observed clustering of clusters. The Ro~D c 
relation for clusters in the low density models, is consistent 
with the scaling relation Ro = 0.4D C for D c < 50ft _1 Mpc. 
On larger scales the cluster correlation length increases more 
slowly than the cluster number density. We also found that 
generally our Ro~D c results were consistent with those of 
Croft & Efstathiou (1994) in spite of the fact that we used 
FOF to identify the clusters in the simulations and Croft 
& Efstathiou (1994) used a very different algorithm. The 
one difference between our results and those of Croft & Ef- 
stathiou (1994) is that in their low Qo run that Ro becomes 
constant for large values of D c whereas we found that Ro 
continues to increase, albeit gently. The analytic Ro~D c re- 
lation based on the PS approach agrees with our numerical 
results. The flattening of the Ro~D c is qualitatively expected 
in CDM models, where the bias of halos depends rather 
weakly on mass (a factor of ~ 2 when the mass changes by 
a factor of 100, see Table |[ while at the same time the num- 
ber density of massive clusters decreases exponentially fast. 
This implies that, as clusters of increasing mass ( i.e. larger 
D c ) are considered, Ro (which is linked to the bias) grows 
slower than requested by the Bachall & West relation. 

Finally, we found that the Ro~D c relation for our low- 
density CDM model is consistent with the results for rich- 
ness R> and R> 1 Abell clusters as well as those for the 
APM clusters, including those for the very rich APM clus- 
ters. The only disagreement is with the single datapoint for 
R> 2 Abell clusters; the correlation length for this sample 
is too large compared with the APM measurements and the 
N-body results. This would suggest some strong systematic 
differences (or selection criteria) between rich APM clusters 
and rich Abell clusters with the same number density. As 
discussed above, it does not seem possible to reconcile CDM 
models with the Bachall & West relation, unless some exotic 
process has boosted the bias of massive clusters compared to 
the dark matter distribution. These considerations hold also 
for flat models with a cosmological constant, as suggested 
by preliminary results by Colberg et al. (1998). 

• Cluster mass function The analytical PS predic- 
tion differs from the N-body mass function both at the 
low and high mass ends (i.e. for both M/M, << 1 and 
M/M, >>1). These discrepancies are more significant in 
the critical SCDM universe. Consequently, the analytic cu- 
mulative cluster temperature function based on the standard 
PS mass function underestimates the abundance of hot clus- 
ters. The discrepancy gets larger with lower as and/or lower 
values of F ( large scale structure studies suggest F ~ 0.25 
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(Maddox et al. 1990)). At z = 1 (assuming a SCDM clus- 
ter normalized model with erg = 0.5) the number of clusters 
with kT> 7 keV is underestimated by almost a factor of ten. 
Claims that rule out critical density models on the basis of 
the detection of massive clusters at high z (Donahue et al. 
1998) should then be taken with some caution. The temper- 
ature function obtained from our simulations evolves more 
slowly than the standard PS temperature function. These 
results however, are not sufficient to bring the cluster tem- 
perature function for flo = 1, CDM models in agreement 
with that observed at z — 0.33 by Henry et al. 1998). 

Our results show there is no evidence to support the 
notion of a universal value for S c in a given cosmology for 
different power spectra. Our results strongly support the 
idea that S c (z) is a function of cosmology, and also of redshift 
and of the mass range under consideration. It is then harder 
to attach a simple physical meaning to the value of S c (z). We 
give a fitting formula for S c (z), valid for groups and clusters 
that improves the agreement between the PS and numerical 
mass functions. It is tempting to interpret the excess of high 
mass objects found in the SCDM simulation as a deviation 
of the gravitational collapse from the idealized spherically 
homogeneous collapse model (Monaco 1995). 

One important, general consideration that can be ex- 
tracted from our analysis is the primary importance of 
comparing results on homogeneous grounds. In particular, 
it is crucial to define cluster samples on the basis of the 
same, physically motivated quantity. If clusters are defined 
in terms of the mass within the virial radius, then both an- 
alytical and numerical results can be directly compared and 
the resulting differences understood. This makes theoretical 
predictions much more robust and easier to compare with 
observational results. 

It appears that collisionless large scale structure simula- 
tions are now in their maturity. Simulations even larger than 
ours are currently being analyzed (see Colberg et al. 1998). 
The statistical results appear to be robust and stable, allow- 
ing a fruitful comparison with observational data, which at 
the same time, are about to experience a manifold growth. 
We are optimistic that a close interplay between theoreti- 
cal and observational results about the internal physics of 
galaxy clusters and their large scale distribution will help to 
improve our knowledge of the physical state of our Universe. 
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